module util_MC
    use rng
    implicit none
    integer, parameter :: naniso=100
    real(8), dimension(0:naniso) :: lt_anisosct, lt_mu
    integer, dimension(100) :: test_mu_s
    real(8) :: max, min
    real(8) :: hsct
    real(8), allocatable, dimension(:) :: pdff, pdfc, cdff, cdfc
    integer, parameter :: MAXFM = 100
    real(8), pointer, dimension(:,:) :: F_form
    real(8), pointer, dimension(:) :: fm_slp
    real(8), pointer, dimension(:,:) :: fm_coef
    integer, parameter :: fphinm=55, fcov1=56, fcov2=57, fcov3=58
    
contains

subroutine init_util_mc
    implicit none
    include "pspec.FI"
    
    allocate(fm_slp(totmesh))
    allocate(fm_coef(2, totmesh))
    
!    open (1234, file = "ktrace.txt", status = 'unknown')
    open (fphinm, file = "phinm.txt", status = 'unknown')
!    open (fcov1, file = "cov1.txt", status = 'unknown')
!    open (fcov2, file = "cov2.txt", status = 'unknown')
!    open (fcov3, file = "cov3.txt", status = 'unknown')
end subroutine

subroutine update_fmshape(psi)
    implicit none
    include "pspec.FI"
    integer :: i
    real(8), pointer, dimension(:) :: psi
    real(8) :: hr
    real(8), dimension(0:totmesh) :: slp
    real(8) :: sr
    hr = 1._8/h(1)
    slp = 0
    do i=1, totmesh-1
        slp(i) = hr*(psi(i+1)-psi(i))
    end do
    do i=1, totmesh
        fm_slp(i) = 0.5_8*(slp(i-1)+slp(i))
        sr = 1._8/(h(i)*psi(i))
        fm_coef(1,i) = fm_slp(i)*sr*0.5
        fm_coef(2,i) = (psi(i)-h(i)*fm_slp(i)*0.5)*sr
    end do
end subroutine

subroutine sample_shape(ntrn, nq, l, nadd, wt)
    use define
    use ifport
    implicit none
    include "pspec.FI"
    integer, intent(in) :: nadd, l
    integer, intent(inout) :: nq
    real(8), intent(in) :: wt
    type(particle), dimension(:), intent(inout) :: ntrn
    integer :: i    
    real(8) :: sp, c, hn, a, b
    sp = 1._8/nadd
    a = fm_coef(1,l)
    b = fm_coef(2,l)
!    if (fm_slp(l)>0._8) then
        do i=1, nadd
            c = -(-0.5_8+i)*sp
            hn = (-b+sqrt(b**2-4*a*c))/(2*a)
            hn = hac(l-1)+hn
!            call addntrn(ntrn, nq, hn, l, 1-2*drand(0), wt, nn_group())
            call addntrn(ntrn, nq, hn, l, 1-2*GetRN(), wt, nn_group())
        end do
!    else
!        do i=1, nadd
!            c = -(-0.5_8+i)*sp
!            hn = (-b+sqrt(b**2-4*a*c))/(2*a)
!            hn = hac(l-1)+hn
!            call addntrn(ntrn, nq, hn, l, 1-2*drand(0), wt, nn_group())
!        end do
!    endif
end subroutine

subroutine sample_flat(ntrn, nq, l, nadd, wt)
    use define
    use ifport
    implicit none
    include "pspec.FI"
    integer, intent(in) :: nadd, l
    integer, intent(inout) :: nq
    real(8), intent(in) :: wt
    type(particle), dimension(:), intent(inout) :: ntrn
    integer :: i    
    real(8) :: sp, c, hn, a, b
    sp = h(l)/nadd
    do i=1, nadd
        hn = hac(l-1) + (-0.5+i)*sp
!        call addntrn(ntrn, nq, hn, l, 1- 2*drand(0), wt, nn_group())
        call addntrn(ntrn, nq, hn, l, 1- 2*GetRN(), wt, nn_group())
    end do
end subroutine

integer function find_mesh(x)
    implicit none
    include "pspec.FI"
    real(8), intent(in) :: x
    integer im
    do im=1, totmesh
        if (x<hac(im)) then
            exit
        end if 
    end do
    find_mesh = im
end function

integer function nn_group()
    use ifport
    implicit none
    include "pspec.FI"
    include "xsec.FI"
    integer g
    real(8) sum, xi
!    xi  = drand(0)
    xi  = GetRN()
    sum = 0.
    do g=1, ng-1
        sum = sum + Xchi(1,g)
        if (xi<sum) then
            exit
        end if 
    end do
    
    nn_group = g
end function 

integer function sample_collision(l, g)
    use ifport
    use define
    implicit none
    include "pspec.FI"
    include "xsec.FI"
    integer, intent(in) :: l, g
    integer tcomp
    real(8) xi, tot, sct
!    xi  = drand(0)
    xi  = GetRN()
    tcomp = complist(l)
    tot = Xtot(tcomp, g)
    sct = Xsct(tcomp, g)
    if (xi>sct/tot) then
        sample_collision = COL_AB
    else 
        sample_collision = COL_SC
    end if 
end function 

!function SCTX_ANI(tcomp, g, mu) result(ret)
!    include "xsec.FI"
!    real(8) :: ret
!    real(8), intent(in) :: mu
!    integer, intent(in) :: tcomp, g
!    
!    ! temporary function
!    ! Sig_n n=1,2,3 = 0.1, n>3 = 0.
!    ret = (Xsct(tcomp, g))
!    
!end function

subroutine init_anisoSct
    integer :: i
    real(8) :: mu
    mu = -1._8
    hsct = 2._8/naniso

    lt_anisosct(0) = 0.
    lt_mu(0) = -1.
    do i=1, naniso
        mu = mu + hsct
        lt_anisosct(i) = anisct_Larsen(mu)
        lt_mu(i) = mu
    end do
    lt_anisosct(naniso) = 1.0
    
    test_mu_s = 0 
    max = -100
    min = 100
end subroutine

function anisct_Larsen(mu) result(val)
    real(8) :: mu
    real(8) :: val
    !val = 0.253183_8 * (1 + mu) * (1.20023 + mu) * (1.52637 - 1.6288*mu + mu**2)
    
    val = 0.255549_8 * (1 + mu) * (1.19651 + mu) * (1.51584 - 1.62508*mu + mu**2)
    !val = 1.16822 *(0.47175+ 0.428*mu - 0.2625*mu**2 + 0.* mu**3 + 0.21875 *mu**4)
end function

function mu_anisotropic result(mu_s)
    use ifport
    real(8) :: xi, mu_s
    integer :: i
!    xi  = drand(0)
    xi  = GetRN()
    do i=0, naniso-1
        if (lt_anisosct(i)>xi) exit
    end do
    mu_s = (xi-lt_anisosct(i-1))/(lt_anisosct(i)-lt_anisosct(i-1))*hsct+ lt_mu(i-1)
!    i = (mu_s+1)*50+1
!    test_mu_s(i) = test_mu_s(i)+1
!    if (mu_s > max) max = mu_s
!    if (mu_s < min) min = mu_s
end function


subroutine sample_direction_anisotropic(mup, mun, dir)
    use define
    use ifport
    implicit none
    real(8), intent(in) :: mup  
    real(8), intent(out) :: mun
    integer, intent(out) :: dir
    real(8) :: xi, phi
    
    mun = mu_anisotropic()
    phi = 2*PI*GetRN() 
    mun = mup*mun - sqrt(1-mup**2)*sqrt(1-mun**2)*cos(phi)
    if (mun>0.) then
        dir = 1
    else
        dir = -1
    end if    
end subroutine

integer function sample_scattering_group(l, g)
    use ifport
    use define
    implicit none
    include "pspec.FI"
    include "xsec.FI"
    integer, intent(in) :: l, g
    integer tcomp, gp
    real(8) sum, xi
!    xi  = drand(0)
    xi  = GetRN()
    sum = 0.
    tcomp = complist(l)
    do gp=1, ng-1
        sum = sum + XsctM(gp, g, tcomp)/Xsct(tcomp, g)
        if (xi<sum) then
            exit
        end if 
    end do
    
    sample_scattering_group = gp
end function 

real(8) function sample_distance(l, g)
    use ifport
    implicit none
    include "xsec.FI"
    include "pspec.FI"
    integer, intent(in) :: l, g
    integer tcomp
    real(8) x_T, xi
!    xi = drand(0)
!    tcomp = complist(l)
!    x_T = Xtot(tcomp, g)
!    sample_distance = -1./x_T * log(xi)
!    sample_distance = -1./Xtot(complist(l), g) * log(drand(0))
    sample_distance = -1./Xtot(complist(l), g) * log(GetRN())
end function

subroutine sample_Nneutron(l, g, w, k, ke, nn, nnw, nwt)
    use ifport
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    integer, intent(in) :: l, g
    real(8), intent(in) :: w, k, ke
    integer, intent(out) :: nn, nnw
    real(8), intent(out) :: nwt
    integer :: tcomp
    real(8) x_T, xi, nuf
    !xi = drand(0)
    xi = GetRN()
    tcomp = complist(l)
    x_T = Xtot(tcomp, g)
    nuf = Xnuf(tcomp, g)
!    nwt = 1._8
!    nn = 1./k * nuf/x_T*w + xi
    nwt = pwadjr(l)
    nn = pwadj(l)/k * nuf/x_T*w + xi
    nnw = 1./ke * nuf/x_T*w + GetRN()
    gk = gk + nuf/x_T*w
!    nfn(l)=nfn(l)+1
!    write(1234, *) gk
end subroutine


subroutine sample_direction(mu, dir)
    use ifport
    implicit none
    real(8), intent(inout) :: mu
    integer, intent(out) :: dir
    real(8) xi
    !xi = drand(0)
    xi = GetRN()
!    if (xi==0. .or. xi==1.) then
!        xi = drand(0)
!    end if 
    mu = 1- 2* xi
!    if (mu.eq.0.)then
!        mu = 1- 2*drand(0)
!    end if
    if (mu>0.) then
        dir = 1
    else
        dir = -1
    end if
end subroutine


subroutine dist2surf(x, l, dir, surf, dist)
    include "pspec.FI"
    integer, intent(in) :: l, dir
    real(8), intent(in) :: x
    integer, intent(out) :: surf
    real(8), intent(out) :: dist
    if (dir>0) then
        dist = hac(l) - x
        surf = l+1
    else
        dist = x - hac(l-1)
        surf = l-1
    end if
end subroutine

subroutine estimateK (nht, k_col, k_trk)
    include "pspec.FI"
    include "xsec.FI"
    include "arrays.FI"
    integer, intent(in)  :: nht
    real(8), intent(out) :: k_col, k_trk
    integer tcomp, i, g
    real(8) nuf, tot
!    k = 0
!    do g=1, ng
!        do i=1, totmesh
!            tcomp = complist(i)
!            nuf = Xnuf(tcomp, g)
!            tot = Xtot(tcomp, g)
!            k = k + nuf*trk_est(i, g)
!        end do
!    end do
!    tot = sum(trk_est)
!    k = k / nht
    gk = gk / nht
    k_col = gk
    k_trk = gktrk/nht
end subroutine

subroutine make_initial_particle(nht, ntrn, bat)
    use define
    use ifport
    implicit none
    include "pspec.FI"
    integer, intent(in) :: nht, bat
    type(particle), dimension(:, :), intent(inout) :: ntrn
    integer i, j, k
    integer :: nq, l
    real(8) :: x, mu
    nq = 0
    do i=1, nht
        !x = drand(0) * sizeprb
        x = GetRN() * sizeprb
        !mu= 1- 2*drand(0)
        mu= 1- 2*GetRN()
        l = find_mesh(ntrn(i, bat)%x)
        call addntrn(ntrn(:,bat), nq, x, l, mu, 1._8, nn_group())
    end do
end subroutine  

subroutine init_pdf
    implicit none
    include "pspec.FI"
    allocate (pdff(0:totmesh))
    allocate (cdff(0:totmesh))
    pdff(0) = 0.
    cdff(0) = 0.
end subroutine

subroutine init_pdf_coarse
    implicit none
    include "pspec.FI"
    allocate (pdfc(0:totm_cmfd))
    allocate (cdfc(0:totm_cmfd))
    pdfc(0) = 0.
    cdfc(0) = 0.
end subroutine


subroutine make_pdf(phi_f)
    implicit none
    include "pspec.FI"
    include "xsec.FI"
    real(8), pointer, dimension(:,:) :: phi_f
    real(8), dimension(totmesh) :: fis
    integer :: i, g
    
    fis = 0.
    do g=1, ng
        do i=1, totmesh
            fis(i) = fis(i) + Xnuf(complist(i), g)*phi_f(i, g)
        end do
    end do
    do i=1, totmesh
        pdff(i) = fis(i)
        cdff(i) = cdff(i-1) + fis(i)
    end do
    pdff = pdff / cdff(totmesh)
    cdff = cdff / cdff(totmesh)
endsubroutine

subroutine make_pdf_coarse(phi_f)
    implicit none
    include "pspec.FI"
    include "xsec.FI"
    real(8), pointer, dimension(:,:) :: phi_f
    real(8), dimension(totmesh) :: fis
    integer :: i, g, ic 
    
    fis = 0.
    pdfc = 0.
    cdfc = 0.
    do g=1, ng
        do i=1, totmesh
            fis(i) = fis(i) + Xnuf(complist(i), g)*phi_f(i, g)
        end do
    end do
    do i=1, totmesh
        ic = fi2co(i)
        pdfc(ic) = pdfc(ic) + fis(i)
    end do
    do i=1, totm_cmfd
        cdfc(i) = cdfc(i-1) + pdfc(i)
    enddo 

    pdfc = pdfc / cdfc(totm_cmfd)
    cdfc = cdfc / cdfc(totm_cmfd)
endsubroutine

subroutine initp_from_dist(nht, ntrn, bat)
    use define
    use ifport
    implicit none
    include "pspec.FI"
    include "arrays.FI"
    include "xsec.FI"
    integer, intent(in) :: nht, bat
    type(particle), dimension(:, :), intent(inout) :: ntrn
    integer i, l, nq
    real(8) ::xi, x, mu
    nq = 0
    do i=1, nht
        !xi = drand(0)
        xi = GetRN()
        do l=1, totmesh-1
            if (xi<cdff(l)) then
                exit
            end if
        end do
        x = GetRN() * h(l) + hac(l-1)
        mu= 1- 2*GetRN()
#ifdef _2NDHRMN
        if (l<totmesh/2) then
            call addntrnm(ntrn(:,bat), nq, x, l, mu, 1._8, -1._8, nn_group())
        else
            call addntrnm(ntrn(:,bat), nq, x, l, mu, 1._8, 1._8, nn_group())
        endif
#else
        call addntrn(ntrn(:,bat), nq, x, l, mu, 1._8, nn_group())
#endif
    end do
end subroutine  

subroutine initp_from_dist_uniform(nht, ntrn, bat, nq)
    use define
    use ifport
    implicit none
    include "pspec.FI"
    include "arrays.FI"
    include "xsec.FI"
    integer, intent(in) :: nht, bat
    integer, intent(out) :: nq
    type(particle), dimension(:, :), intent(inout) :: ntrn
    integer i, j, intrn
    real(8) :: wtloc, delh, w, x, mu
    integer :: srcloc
    intrn = 1
    nq = 0
    do i=1, totmesh
        wtloc = pdff(i) * nht
        srcloc = nint(wtloc)
        delh = h(i)/srcloc
        if (srcloc .eq. 0) cycle
        if (swfmshape) then
            call sample_shape(ntrn(:,bat), nq, i, srcloc, wtloc/srcloc)
        else
            call sample_flat(ntrn(:,bat), nq, i, srcloc, wtloc/srcloc)
        endif
!        do j=1, srcloc
!!            ntrn(intrn, bat)%x  = hac(i-1) + (-0.5 + j) * delh
!!            ntrn(intrn, bat)%mu = 1- 2*drand(0)
!!            ntrn(intrn, bat)%w  = wtloc/srcloc
!!            ntrn(intrn, bat)%l  = i
!!            ntrn(intrn, bat)%g  = nn_group()
!!            intrn = intrn+1
!            x = hac(i-1) + (-0.5 + j) * delh
!            mu= 1- 2*drand(0)
!            w  = wtloc/srcloc
!            call addntrn(ntrn(:,bat), nq, x, i, mu, w, nn_group())
!        end do
    end do
!    nq = intrn-1
end subroutine


subroutine initp_from_dist_uniform_coarse(nht, ntrn, bat, nq)
    use define
    use ifport
    implicit none
    include "pspec.FI"
    include "arrays.FI"
    include "xsec.FI"
    integer, intent(in) :: nht, bat
    integer, intent(out) :: nq
    type(particle), dimension(:, :), intent(inout) :: ntrn
    integer i, j, l, intrn
    real(8) :: wtloc, delh, x, mu, w
    integer :: srcloc, nfine
    nfine = totmesh/totm_cmfd
    intrn = 1
    nq = 0 
    do i=1, totm_cmfd
        wtloc = pdfc(i) * nht
        srcloc = nint(wtloc)
        delh = h_cmfd(i)/srcloc
        l=1
        do j=1, srcloc
!            ntrn(intrn, bat)%x  = hac_cmfd(i-1) + (-0.5 + j) * delh
!            if (ntrn(intrn, bat)%x>hac(nfine* (i-1) + l)) l=l+1
!            ntrn(intrn, bat)%mu = 1- 2*drand(0)
!            ntrn(intrn, bat)%w  = wtloc/srcloc
!            ntrn(intrn, bat)%l  = nfine* (i-1) + l
!            ntrn(intrn, bat)%g  = nn_group()
!            intrn = intrn+1
            x = hac_cmfd(i-1) + (-0.5 + j) * delh
            if (x>hac(nfine* (i-1) + l)) l=l+1
!            mu= 1- 2*drand(0)
            mu= 1- 2*GetRN()
            w  = wtloc/srcloc
            call addntrn(ntrn(:,bat), nq, x, nfine* (i-1) + l, mu, w, nn_group())
       end do
    end do
!    nq = intrn-1
end subroutine



subroutine initp_from_coarsedist(nht, ntrn, bat)
    use util_CMFD
    use define
    use ifport
    implicit none
    include "pspec.FI"
    include "arrays.FI"
    include "xsec.FI"
    include "hmgz.FI"

    integer, intent(in) :: nht, bat
    type(particle), dimension(:, :), intent(inout) :: ntrn
    real(8), dimension(totm_cmfd) :: fis
    real(8), dimension(0:totm_cmfd) :: pdf
    integer i, g, l, j, nq
    real(8) ::xi, x, mu, w
    
    fis = 0.
    pdf(0) = 0.
    nq = 0
    do g=1, ng
        do i=1, totm_cmfd
            fis(i) = fis(i) + hNuf(i, g)*hphi(i, g)
        end do
    end do
    do i=1, totm_cmfd
        pdf(i) = pdf(i-1) + fis(i)
    end do
    pdf = pdf / pdf(totm_cmfd)
    
    do i=1, nht
!        xi = drand(0)
        xi = GetRN()
        do l=1, totm_cmfd-1
            if (xi<pdf(l)) then
                exit
            end if
        end do
!        ntrn(i, bat)%x  = drand(0) * h_cmfd(l) + hac_cmfd(l-1)
!        ntrn(i, bat)%mu = 1- 2*drand(0)
!        ntrn(i, bat)%w  = 1.
!        ntrn(i, bat)%l  = l
!        ntrn(i, bat)%g  = nn_group()
!        x = drand(0) * h_cmfd(l) + hac_cmfd(l-1)
        x = GetRN() * h_cmfd(l) + hac_cmfd(l-1)
!        mu= 1- 2*drand(0)
        mu= 1- 2*GetRN()
        w  = 1._8
        call addntrn(ntrn(:,bat), nq, x, l, mu, w, nn_group())
    end do
end subroutine  

subroutine add_particle(nq, ntrn, bat, nadd)
    use util_CMFD
    use define
    use ifport
    implicit none
    include "pspec.FI"
    include "arrays.FI"
    include "xsec.FI"
    integer, intent(inout) :: nq
    integer, intent(in) :: bat, nadd
    type(particle), dimension(:, :), intent(inout) :: ntrn
    real(8) :: x
    real(8) :: a, b, del
    integer :: i,l
    integer :: diagram(totmesh)
    a = 5.; b= 10.
    del = b-a
    diagram = 0
    do i=1, nadd
!        ntrn(nq+i, bat)%x  = a + del/PI*(acos(1-2*drand(0)))
        ntrn(nq+i, bat)%x  = a + del/PI*(acos(1-2*GetRN()))
!        ntrn(nq+i, bat)%mu = 1- 2*drand(0)
        ntrn(nq+i, bat)%mu = 1- 2*GetRN()
        ntrn(nq+i, bat)%w  = 1.
        do l=1, totmesh
            if (ntrn(nq+i, bat)%x<hac(l)) exit
        end do
        diagram(l-1) = diagram(l-1)+1
        ntrn(nq+i, bat)%l  = l-1
        ntrn(nq+i, bat)%g  = nn_group()
    end do
    
    nq = nq + nadd
    
end subroutine


subroutine leak2netj(nht)  
    use allocs  
    use define
    implicit none
    include "arrays.FI"
    include "accm.FI"
    include "pspec.FI"
    include "xsec.FI"
    integer, intent(in) :: nht
    integer :: i, g
    real(8) :: netj
    real(8), dimension(totmesh*2) :: tempJ
    jf = jf/nht
    pjf = jf
    do g=1, ng
        tempJ = jf(:, g)
        do i=2, totmesh*2 -1, 2
            netj = (tempJ(i)-tempJ(i+1))
            jf(i,g) = netj
            jf(i+1,g) = -netj
        end do
    end do
    
    jfAC = jfAC + jf
    pjAC = pjAC + pjf
end subroutine

subroutine evaluate_flux (nht,ifst)
    implicit none
    include "arrays.FI"
    include "accm.FI"
    include "pspec.FI"
    include "xsec.FI"
    include "stsol.FI"
    include "CMFD_LARSEN.FI"
    integer, intent(in) :: nht
    integer :: g, i
    real(8) :: r, l2
    logical :: ifst
    integer :: nfp
    real(8) :: vpsi, cyck
    do g = 1, ng
        do i=1, totmesh
            phi_f(i, g) = trk_est(i, g)/nht/h(i)
            phi_f2(i,g) = phi_f2(i,g)/nht/h(i)
#ifdef _2NDHRMN 
            phi_f2nd(i,g) = trk_est2(i,g)/nht/h(i)
#endif
        end do
    end do
    
#ifdef CMFDLARSEN
    do g=1, ng
        do i=1, totmesh
            phi1La(i,g) = phi1La(i,g)/nht/h(i)
        enddo
    enddo
    phi2Ls = phi2Ls / nht
#endif    
    phifAC = phifAC + phi_f
    phif2AC = phif2AC + phi_f2
    psi_f = 0.
    psi_f2 = 0.
    do g=1, ng
        do i=1, totmesh
            psi_f(i) = psi_f(i) + Xnuf(complist(i), g) * phi_f(i, g) * h(i)
            ! for 2nd harmonics            
#ifdef _2NDHRMN
            psi_f2(i) = psi_f2(i) + Xnuf(complist(i), g) * trk_est2(i,g)/nht 
#endif
        enddo
    enddo
    
    nfp = 0
    do i=1, totmesh
        if (psi_f(i)>0._8) nfp=nfp+1
    enddo
    
    
    r = 1._8/sum(psi_f)*sizeprb
    phi_fnm = phi_f*r
    
    if (ifst) then
    cyck=gk/nht
    do g=1, ng
        stmc(0,g)=stmc(0,g)+cyck
        stsqmc(0,g)=stsqmc(0,g)+cyck**2
        do i=1, totmesh
            stmc(i,g) = stmc(i,g) + phi_fnm(i,g)
            stsqmc(i,g) = stsqmc(i,g) + (phi_fnm(i,g))**2
        enddo
    enddo
    
    r=r/h(1)*nfp/totmesh
    stmcf(0)=stmcf(0)+cyck
    stsqmcf(0)=stsqmcf(0)+cyck**2
    do i=1, totmesh
        vpsi=psi_f(i)*r
        stmcf(i)=stmcf(i)+vpsi
        stsqmcf(i)=stsqmcf(i)+vpsi**2
    enddo
        
    endif


    r=1._8/sum(stmc)
    err_ref=0.
    do g=1, 1
        do i=1, totmesh
            err_ref=err_ref+(refflux(i,g)-stmc(i, g)*r)**2
        enddo
    enddo
    err_ref=sqrt(err_ref*totmesh)
    
    
#ifndef _NOLOG
    write(fphinm, '(2000(e11.5,1x))') (phi_fnm(i,1),i=1,totmesh)
#endif
    
    do i=1, totmesh
        psifAC(i) = psifAC(i) + psi_f(i)
    enddo

    do g=1, ng
        do i=1, totmesh+1
            phi_s (i,g)= phi_s (i,g)/nht
            mmt1st(i,g)= mmt1st(i,g)/nht
            mmt2nd(i,g)= mmt2nd(i,g)/nht
         enddo
         do i=1, totm_cmfd+1
            phi_sa(i,g)= phi_sa(i,g)/nht
            mmtAbsr(i,g)= mmtAbsr(i,g)/nht
            mmtNuF (i,g)= mmtNuF (i,g)/nht
        enddo
    enddo
    phisAC   = phisAC   + phi_s
    phisaAC  = phisaAC  + phi_sa
    mmt1stAC = mmt1stAC + mmt1st
    mmt2ndAC = mmt2ndAC + mmt2nd
    mmtAbsrAC = mmtAbsrAC + mmtAbsr
    mmtNuFAC  = mmtNuFAC  + mmtNuF
!    if (PDFREF) then
!        l2=0
!        do g=1, ng
!            do i=1, totmesh
!                l2 = l2 + (refphi(i,g)-phi_f(i,g))**2
!            end do
!        end do
!        l2 = sqrt(l2)/l2ref_f
!        gtally(1) = l2
!    end if
    if (PDFREF) call L2ERR
end subroutine

subroutine L2ERR
    include "pspec.FI"
    include "xsec.FI"
    include "arrays.FI"
    include "accm.FI"
    real(8) :: nf, l2
    integer :: i
    l2 = 0 
    nf = 1._8/sum(psifac)
    do i=1, totmesh
        l2 = l2+(refpsi(i)-psifac(i)*nf)**2
    end do
    l2 = sqrt(l2)/l2ref_f
    gtally(1) = l2
end subroutine

subroutine init_weight(nqn, nwsum, ntrn)
    use define
    implicit none
    include "arrays.FI"
    include "pspec.FI"
    include "xsec.FI"
    integer, intent(in) :: nqn, nwsum
    type(particle), dimension(:), intent(inout) :: ntrn
    integer ::i, l, g
    real(8) w, wsum, ratio
    
    do i=1, nqn
        wsum = wsum+ntrn(i)%w
    end do
    
    ratio=(1._8*nwsum/wsum) !*(1._8*nwsum/nqn)
    do i=1, nqn
        ntrn(i)%w = ntrn(i)%w*ratio
    end do
    
    popctrl=1._8*nwsum/nqn
    
end subroutine

subroutine feedback(nqn, nwsum, ntrn)
    use define
    implicit none
    include "arrays.FI"
    include "pspec.FI"
    include "xsec.FI"
    integer, intent(in) :: nqn, nwsum
    type(particle), dimension(:), intent(inout) :: ntrn
    integer ::i, ic, l, g
    real(8) w, wsum, ratio
    real(8) :: nsrc(totm_cmfd)
    real(8) :: wtadj(totm_cmfd)
    nsrc = 0
    do i=1, nqn
        l = ntrn(i)%l
        ic = fi2co(l)
        nsrc(ic) = nsrc(ic)+ntrn(i)%w
    end do
    wtadj = 1
    
!    if (PDFREF) then
!        do i=1, totm_cmfd
!            if (nsrc(i)>0) wtadj(i) = refPDF(i) * nwsum / nsrc(i)
!        end do
!    else
        do i=1, totm_cmfd
            if (nsrc(i)>0) wtadj(i) = srcPDF(i) * nwsum / nsrc(i)
        end do
!    endif
    
    do i=1, nqn
        l = ntrn(i)%l
        ic = fi2co(l)
        ntrn(i)%w = ntrn(i)%w * wtadj(ic)
    end do
end subroutine

subroutine feedbackRSMP(nqn, nwsum, ntrn)
    use define
    implicit none
    include "arrays.FI"
    include "pspec.FI"
    include "xsec.FI"
    integer, intent(inout) :: nqn
    integer, intent(in) :: nwsum
    type(particle), dimension(:), intent(inout) :: ntrn
    integer ::i, j, l, g
    real(8) :: w, wsum, xi, x
    real(8) :: srcCDF(0:totm_cmfd)
    real(8) :: check(totm_cmfd)
    integer :: nr
    
    srcCDF = 0
!    if (PDFREF) then
!        do i=1, totm_cmfd
!            srcCDF(i) = srcCDF(i-1) + refPDF(i)
!        end do
!    else
        do i=1, totm_cmfd
            srcCDF(i) = srcCDF(i-1) + srcPDF(i)
        end do
!    endif        
    nqn = 0
    check = 0
    do i=1, nwsum
        xi = getRN()
        do j=1, totm_cmfd-1
            if (xi< srcCDF(j)) exit
        end do
        x = (xi-srcCDF(j-1))/(srcCDF(j)-srcCDF(j-1))*h_cmfd(j)+hac_cmfd(j-1)
        nr = totmesh/totm_cmfd
        if (nr>1) then
            do l=(j-1)*nr, j*nr-1
                if (x<hac(l)) exit
            enddo
        else 
            l = j
        endif
        call addntrn(ntrn, nqn, x, l, 1-2*getRN(), 1._8, nn_group())    
        check(j) = check(j)+1
    end do
end subroutine

subroutine feedback_FMC(nqn, nwsum, ntrn)
    use define
    implicit none
    include "arrays.FI"
    include "pspec.FI"
    include "xsec.FI"
    integer, intent(in) :: nqn, nwsum
    type(particle), dimension(:), intent(inout) :: ntrn
    integer ::i, ic, l, g
    real(8) w, wsum, ratio, hh
    integer :: nsrc(totmesh+1)
    real(8) :: wtadj(totmesh+1), sumt
    nsrc = 0
    do i=1, nqn
        l = ntrn(i)%l
        hh = hac(l)-ntrn(i)%x
        if (hh> h(l)/2) then
            nsrc(l) = nsrc(l)+1
        else
            nsrc(l+1) = nsrc(l+1)+1
        endif 
    end do
    do i=1, totmesh+1
        if (nsrc(i)>0) wtadj(i) = adj_fmc(i) * nwsum / nsrc(i)
    end do
    
!    sumt  = 0;
    do i=1, nqn
        l = ntrn(i)%l
        hh = hac(l)-ntrn(i)%x
        if (hh> h(l)/2) then
            ntrn(i)%w = ntrn(i)%w * wtadj(l)
        else
            ntrn(i)%w = ntrn(i)%w * wtadj(l+1)
        endif 
!        sumt = sumt + ntrn(i)%w
    end do
end subroutine

subroutine formfunction(phi_f)
    implicit none
    include "pspec.FI"
    include "xsec.FI"
    real(8), dimension(totmesh, ng), intent(in) :: phi_f
    real(8) :: fis(totmesh)
    integer :: nfine 
    integer :: i, j, fi, g
    real(8) :: lcsum
    nfine = totmesh/totm_cmfd
    allocate(F_form(nfine, totm_cmfd))
    fis = 0
    do g=1, ng
        do i=1, totmesh
            fis(i) = fis(i) + Xnuf(complist(i), g)*phi_f(i, g)
        end do
    end do
    fi = 1
    do i=1, totm_cmfd
        lcsum =0;
        do j=1, nfine
            F_form(j, i)= fis(fi)
            lcsum = lcsum + fis(fi)
            fi = fi+1  
        end do
        if (lcsum .ne. 0.) F_form(1:nfine,i) = F_form(1:nfine,i)/lcsum
    end do
end subroutine

subroutine update_pdf_coarse(psi_c)
    implicit none
    include "pspec.FI"
    real(8), dimension(totmesh) :: psi_c
    integer :: i
    
    do i=1, totm_cmfd
        pdfc(i) = psi_c(i)
        cdfc(i) = cdfc(i-1) + pdfc(i)
    enddo 

    pdfc = pdfc / cdfc(totm_cmfd)
    cdfc = cdfc / cdfc(totm_cmfd)    
    
end subroutine

subroutine initp_from_shapeF(nht, ntrn, bat, nq)
    use define
    use ifport
    implicit none
    include "arrays.FI"
    include "pspec.FI"
    integer, intent(in) :: nht, bat
    integer, intent(out) :: nq
    type(particle), dimension(:, :), intent(inout) :: ntrn
    real(8) :: wtloc, wtloc2 
    integer i, j, k, l, fi, nfine, intrn, srcloc
    real(8) :: delh, x, mu, w
    nfine = totmesh/totm_cmfd
    fi = 1
    intrn = 1
    nq = 0
    do i=1, totm_cmfd
        wtloc = pdfc(i) * nht
        do j=1, nfine
            wtloc2 = wtloc * F_form(j,i)
            srcloc = nint(wtloc2)
            delh = h(i)/srcloc
            if (swfmshape) then
                call sample_shape(ntrn(:,bat), nq, fi, srcloc, wtloc2/srcloc)
            else
                call sample_flat(ntrn(:,bat), nq, fi, srcloc, wtloc2/srcloc)
            endif
!            do k=1, srcloc
!!                ntrn(intrn, bat)%x  = hac(fi-1) + (-0.5 + k) * delh
!!                ntrn(intrn, bat)%mu = 1- 2*drand(0)
!!                ntrn(intrn, bat)%w  = wtloc2/srcloc
!!                ntrn(intrn, bat)%l  = fi
!!                ntrn(intrn, bat)%g  = nn_group()
!!                intrn = intrn+1
!                x = hac(fi-1) + (-0.5 + k) * delh
!                mu= 1- 2*drand(0)
!                w  = wtloc2/srcloc
!                call addntrn(ntrn(:,bat), nq, x, l, mu, w, nn_group())
!            end do
            fi = fi+1
        end do
    end do
!    nq = intrn-1
end subroutine
!subroutine feedback_shape(wtsum, np, ntrn)
!    use define
!    implicit none
!    include "arrays.FI"
!    include "pspec.FI"
!    include "xsec.FI"
!    integer, intent(in) :: wtsum
!    integer, intent(inout) :: np
!    type(particle), dimension(:), intent(inout) :: ntrn
!    integer ::i, j, l, g, nfine
!    real(8) w, wsum, ratio
!    real(8) :: pdfc(0:totm_cmfd)
!    real(8) :: pdfl(MAXFM)
!    nfine = totmesh/totm_cmfd
!    if (nfine>MAXFM) then
!        write(*) "ERR : Too many fine-meshes in a coarse-mesh!"   
!        exit
!    endif
!    ! make pdf for coarse meshes
!    pdfc = 0
!    pdfc(1:totm_cmfd) = psi_c(1:totm_cmfd)
!    
!    do i=1. totm_cmfd
!        do j=1, nfine
!        end do
!    end do
!    do i=1, nht
!        w = ntrn(i)%w
!        l = ntrn(i)%l
!        g = ntrn(i)%g
!        ntrn(i)%w = w * adjp(l)
!    end do
!
!    wsum = 0.    
!    do i=1, nht 
!        wsum = wsum + ntrn(i)%w
!    end do
!    ratio = 1._8/wsum*np
!!    ratio = 1./ratio
!    do i=1, nht
!        ntrn(i)%w = ntrn(i)%w*ratio
!    end do
!    
!end subroutine


subroutine adjustWeight(wtsum, nHistory, ntrn, nptc, wtm)
!   This routine won't work well when implicit captuer scheme is applied. 
!   Because it adjust weighting just before generating new source neutrons
!   That means just before changing # of neutrons, each neutron has adjusted weight.
!   But it dosen't matter for the results of tallies, since the wtsum is preserved. 
    use define
    implicit none
    include "arrays.FI"
    include "pspec.FI"
    include "xsec.FI"
    real(8), intent(inout) :: wtsum
    integer, intent(inout) :: nHistory
    integer, intent(in) :: nptc
    real(8), intent(in) :: wtm
    type(particle), dimension(:), intent(inout) :: ntrn
    integer ::i
    
    wtsum = wtsum*wtm
    nHistory = nHistory*wtm
    
    do i=1, nptc
        ntrn(i)%w = ntrn(i)%w * wtm
    end do
end subroutine

subroutine strcat(str1, str2)
    character*200:: str1,str2

end subroutine

function ShannonEntropy
    implicit none
    include "arrays.FI"
    include "pspec.FI"
    include "xsec.FI"
    real(8) :: ShannonEntropy, pj, sumpsi
    integer :: i, g, fismesh
	ShannonEntropy = 0.
	fismesh =0
    sumpsi = sum(psi_f)
!	do i=1, totmesh
!	    if (psi_f(i)>0.) fismesh = fismesh+1
!	end do
	
	do i=1, totmesh
!	    pj = psi_f(i)/fismesh
	    pj = psi_f(i)/sumpsi
	    if (pj>0.) ShannonEntropy = ShannonEntropy -pj*log(pj)
	end do
	ShannonEntropy = ShannonEntropy/ log(2.)
end function


subroutine addntrn(ntrn, nq, x, l, mu, w, g)
    use define
    include "pspec.FI"
    type(particle), dimension(:), intent(inout) :: ntrn
    integer, intent(inout) :: nq
    integer, intent(in) :: l, g
    real(8), intent(in) :: x, mu, w
    real(8) :: mum
    integer :: i, j
    mum = mu
    do j=1, swsplitting
        nq = nq +1                                                   
        ntrn(nq)%x = x                                          
        ntrn(nq)%w = w                                         
        ntrn(nq)%l = l                                          
        ntrn(nq)%g = g
#ifdef ISORAND    
        ntrn(nq)%mu = 1-2* getRN()
#else
        ntrn(nq)%mu = mum       
        mum = mum + 2._8/swsplitting
        if (mum>1._8) mum = mum-2
#endif
    end do
end subroutine

subroutine addntrnm(ntrn, nq, x, l, mu, w, wm, g)
    use define
    include "pspec.FI"
    type(particle), dimension(:), intent(inout) :: ntrn
    integer, intent(inout) :: nq
    integer, intent(in) :: l, g
    real(8), intent(in) :: x, mu, w, wm
    real(8) :: mum
    integer :: i, j
    mum = mu
    do j=1, swsplitting
        nq = nq +1                                                   
        ntrn(nq)%x = x                                          
        ntrn(nq)%w = w                                         
        ntrn(nq)%wm = wm                                         
        ntrn(nq)%l = l                                          
        ntrn(nq)%g = g
#ifdef ISORAND    
        ntrn(nq)%mu = 1-2* getRN()
#else
        ntrn(nq)%mu = mum       
        mum = mum + 2._8/swsplitting
        if (mum>1._8) mum = mum-2
#endif
    end do
end subroutine

subroutine switch_cycle(bat)
    integer :: bat(3)
    integer :: temp
    temp = bat(1)
    bat(1) = bat(2)
    bat(2) = temp
end subroutine
subroutine switch_wielandt(bat)
    integer :: bat(3)
    integer :: temp
    temp = bat(1)
    bat(1) = bat(3)
    bat(3) = temp    
end subroutine
function shiftk(iter, delk) result (ret)
    real(8) :: SHIFT_GUIDe(18)
    data SHIFT_GUIDE /1.0_8, 1.0_8, 0.75_8, 0.50_8, 0.35_8, 0.25_8, &
                    0.20_8, 0.15_8, 0.125_8, 0.100_8, 0.075_8, 0.050_8, &
                    0.035_8, 0.025_8, 0.020_8, 0.015_8, 0.0125_8, 0.0100_8/
    integer, intent(in) :: iter
    real(8), intent(in) :: delk
    real(8) :: ret
    ret = delk
    if (iter .le. 18) ret = SHIFT_GUIDE(iter) 
end function

!function switch(swt) result(ret)
!    integer :: LT_SWT(2)
!    data LT_SWT /2, 1/
!    integer, intent(in) :: swt
!    integer :: ret
!    ret = LT_SWT(swt)
!end function

subroutine evaluate_corcoef(batch, val1, val2, val3)
    use ifport
    use util_STT
    include "pspec.FI"
    integer :: i,j,k
    integer, intent(in) :: batch
    integer :: st, ed
    integer, parameter :: n=100
    real(8), intent(in) :: val1, val2, val3
    real(8) :: mean(3), Cr(3), sum2(3), rk, e1, e2, s2(3)
    real(8) :: covlag, acrc(n-1,3), stds, stdt
    gtally_arr(1, batch) = val1
    gtally_arr(2, batch) = val2
!    gtally_arr(3, batch) = val3
    gtally_arr(3, batch) = drand(0) !val3

    st = batch-n+1
    ed = st + n - 1
    if (st<1) return
    sum2 = 0
    do i=1, n-2
        Cr=0
        do k=1, 3
            stds = stdevp(gtally_arr(k,:),st,ed-i)
            stdt = stdevp(gtally_arr(k,:),st+i,ed)
            covlag = covarlag(gtally_arr(k,:), st, ed, i)
            acrc(i,k) = covlag/(stds*stdt)
            sum2(k)  = sum2(k) + (n-i)* covlag
        end do 
    end do
#ifdef _COEFF
#ifndef _NOLOG
    write(fcov1,   '(2000(e13.5,1x))') (acrc(i,1),i=1,n-2)    
    write(fcov2,   '(2000(e11.5,1x))') (acrc(i,2),i=1,n-2)    
    write(fcov3,   '(2000(e11.5,1x))') (acrc(i,3),i=1,n-2)    
#endif
#endif
    do k=1, 3
        s2(k) = var(gtally_arr(k,:),st, ed)
    end do
    s2 = s2/n
    gtally = 0
    forall(k=1:3) gtally(k) = sqrt(abs(2._8/n/(n-1)*sum2(k) / s2(k)))
end subroutine

subroutine find2ndHarmonic
    include "pspec.FI"
    include "arrays.FI"
    
    
end subroutine



subroutine write_result
    use statistics, only : avgk_mc, stdk_mc
    include "arrays.FI"
    include "xsec.FI"
    include "pspec.FI"
    include "stsol.FI"
    integer,parameter :: favg=901, fstdev=902
    integer :: g, j, i
    real(8) :: rnr
    character(256) :: sfile1, sfile2
    real(8), dimension(0:totmesh+1) :: stdevbuf
    rnr = 1._8

    open(favg,file ="avgmc.txt",status='unknown')
    open(fstdev,file ="stdmc.txt",status='unknown')
    do g=1, ng
        do i=0, totmesh
            stdevbuf(i) = rnr*sqrt((stsqmc(i,g)-stmc(i,g)**2/active)/active/(active-1))
            stmc(i,g) = stmc(i,g)*rnr/active
        enddo
        write(favg, '(2000(es15.8,1x))') avgk_mc, stmc(1:totmesh,g)
        write(fstdev, '(2000(es15.8,1x))') stdk_mc, stdevbuf(1:totmesh)
    enddo
    close(favg)
    close(fstdev)

    open(favg,file ="avgmcf.txt",status='unknown')
    open(fstdev,file ="stdmcf.txt",status='unknown')
    do i=0, totmesh
        stdevbuf(i) = rnr*sqrt((stsqmcf(i)-stmcf(i)**2/active)/active/(active-1))
        stmcf(i) = stmcf(i)*rnr/active
    enddo
    write(favg, '(2000(es15.8,1x))') avgk_mc, stmcf(1:totmesh)
    write(fstdev, '(2000(es15.8,1x))') stdk_mc, stdevbuf(1:totmesh)
    close(favg)
    close(fstdev)

    if (NMVWND .eq. 1) then
        open(favg,file ="avgcmfd.txt",status='unknown')
        open(fstdev,file ="stdcmfd.txt",status='unknown')
        do g=1, ng
            do i=0, totm_cmfd
                stdevbuf(i) = rnr*sqrt((stsqcmfd(i,g)-stcmfd(i,g)**2/active)/active/(active-1))
                stcmfd(i,g) = stcmfd(i,g)*rnr/active
            enddo
            write(favg, '(2000(e15.8,1x))') stcmfd(0:totm_cmfd, g)
            write(fstdev, '(2000(e15.8,1x))') stdevbuf(0:totm_cmfd)
        enddo
        close(favg)
        close(fstdev)
        
        open(favg,file ="avgcmfdf.txt",status='unknown')
        open(fstdev,file ="stdcmfdf.txt",status='unknown')
        do i=0, totm_cmfd
            stdevbuf(i) = rnr*sqrt((stsqcmfdf(i)-stcmfdf(i)**2/active)/active/(active-1))
            stcmfdf(i) = stcmfdf(i)*rnr/active
        enddo
        write(favg, '(2000(e15.8,1x))') stcmfdf(0:totm_cmfd)
        write(fstdev, '(2000(e15.8,1x))') stdevbuf(0:totm_cmfd)
        close(favg)
        close(fstdev)
        

!        open(favg,file ="avgfmc.txt",status='unknown')
!        open(fstdev,file ="stdfmc.txt",status='unknown')
!        do i=0, totm_cmfd+1
!            stdevbuf(i) = rnr*sqrt((stsqfmc(i)-stfmc(i)**2/active)/active/(active-1))
!            stfmc(i) = stfmc(i)*rnr/active
!        enddo
!        write(favg, '(2000(e15.8,1x))')  stfmc(0:totm_cmfd+1)
!        write(fstdev, '(2000(e15.8,1x))') stdevbuf(0:totm_cmfd+1)
!        close(favg)
!        close(fstdev)
    else 
        open(favg,file ="avgcmfd.txt",status='unknown')
        do g=1, ng
            do i=0, totm_cmfd
                stcmfd(i,g) = stcmfd(i,g)
            enddo
            write(favg, '(2000(e15.8,1x))') stcmfd(0:totm_cmfd, g)
        enddo
        close(favg)
        open(favg,file ="avgcmfdf.txt",status='unknown')
        do i=0, totm_cmfd
            stcmfdf(i) = stcmfdf(i)
        enddo
        write(favg, '(2000(e15.8,1x))') stcmfdf(0:totm_cmfd)
        close(favg)
    endif   
    
    open(favg,file ="nfn.txt",status='unknown')
    write(favg, '(2000(i8,1x))') nfn
    close(favg)
    

end subroutine

end module

